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Aims. We investigate the dynamics of a circumbinary disc that responds to the loss of mass and to the recoil velocity of the black hole 
produced by the merger of a binary system of supermassive black holes. 

Methods. We perform the first two-dimensional general relativistic hydrodynamics simulations of extended non-Keplerian discs and 
employ a new technique to construct a "shock detector", thus determining the precise location of the shocks produced in the accreting 
disc by the recoiling black hole. In this way we can study how the properties of the system, such as the spin, mass and recoil velocity 
of the black hole, affect the mass accretion rate and are imprinted on the electromagnetic emission from these sources. 
Results. We argue that the estimates of the bremsstrahlung luminosity computed without properly taking into account the radiation 
transfer yield cooling times that are unrealistically short. At the same time we show, through an approximation based on the relativistic 
isothermal evolution, that the luminosity produced can reach a peak value above L = 10 43 erg/s at about ~ 30 d after the merger of 
a binary with total mass M a 10 6 M G and persist for several days at values which are a factor of a few smaller. If confirmed by 
more sophisticated calculations such a signal could indeed lead to an electromagnetic counterpart of the merger of binary black-hole 
system. 
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1. INTRODUCTION 

Despite the fundamental role they play in gravitational-wave as- 
tronomy, no undisputed observational evidence of the existence 
of supermassive binary black holes (SMBBHs) systems has been 
found yet. However, circumstantial evidence does exist for a 
number of potential candidates. This is the case, for instance, 
of the radio galaxy 0402+379, which shows a projected separa- 
tion between th e two black holes of 7 .3 pc and a total mass of 
~ 1.5 x 10 8 M o (Rodrigu ez et al.|2006) . Similarly, the ultralumi- 
nous infrared galaxy NGC6240 shows two optical nuclei and is 
thought to be in an early merger phase (Komo ssa et al.||2"003) . 
Finally, potential candidates have been suggested in few other 
cases where the two galaxies are more widely separated, like in 
the pair IC694/NGC3690 hosting two active nuclei as revealed 
by the presence of two distinct Ka lines in their X-ray spectra 
HBallo et al.|2004) . 

In addition, there is a large family of even more uncertain 
SMBBH candidates, that are spatially unresolved and whose ul- 
timate nature is, of course, a matter of strong debate. They in- 
clude the class of X-shaped radio galaxies, in which the observed 
changes in the orientation of the black hole spin axis could be 
due to an ongoing merger with a second black hole (Gopal- 
|Krishna et al.|2 003); or the class of double-double radio galax- 
ies presenting a pair of double-lobed radio structures that could 
be the remnants of a SMBBH merger event ( Liu et aL|2003| l; 
or, finally, the class of sources showing periodicities in the light 
curves, like in BL Lac Object OJ287 ( |Komo ssa 2006 ). Quite re- 
cently, also the quasar SDSSJ0927 (at a redshift z ~ 0.7) has 
been identified as a promising SMBBH candidate, with a mass 



for the primary black hole of M « 2 x 10 9 M Q and a semi-major 
axis of 0.34 pc ( jDotti et alT2009 1. 



A strong motivation for studying the merger of super- 
massive binary black hole systems comes from the fact that 
their gravitational signal will be detected by the planned Laser 
Interferometric Space Antenna (LISA), whose optimum sensi- 
tivity is placed in the range (10~ 4 -r 0.1) Hz. Considerable at- 
tention has therefore been recently attracted by the possibility 
of detecting also the electromagnetic (EM) counterpart of these 
events through the emission coming from the circumbinary ac- 
cretion disc that is expected to form when the binary is still 
widely separated. Such a detection will not only act as a con- 
firmation of the gravitational wave (GW) detection, but it will 
also provide a new tool for testing a number of fundamental as- 
trophysical issues (Hai marTet al.|20 09). More specifically, it will 
offer the possibility of testing models of galaxy mergers and ac- 
cretion discs perturbation, probing basic aspects of gravitational 
physics, and allowing for the measurement of the Eddington ra- 
tio and for the determination of cosmological parameter once 
the redshift is known (Phinney 2009). In spite of the presence of 
the disturbing effects due to weak-lensing errors, |Kocsis et al.| 
(2006 1 have computed the average number of quasars in the 
three-dimensional LISA error volume and have shown that for 
mergers with masses in the range ~ 4 x (10 5 -r 1O 7 )M at red- 
shift z ~ 1, the error box may cont ain 1 quasar with a l uminosity 
L B ~ (10 10 h- 10 n )L G (see also |Kocsis et~al|2008| [Kocsis & 
|Loeb|2008) . 

As a product of this increased interest, a number of studies 
have been recently carried out to investigate the properties of 
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these EM counterparts either during the stages that precede the 
merger, or in those following it. As an example, recent work has 
considered the interaction between the binary and a dense gas 
cloud ([Armitage & Natarajan|2002"t|van Meter et al.|2009}|Bode 
et1QT|2009| |Farris et al.||2009| |Lodato et al.||2009| [Chang et aT 



2009 ) even though astrophysical considerations seem to suggest 



that during the very final stages of the merger the SMBBH will 
inspiral in a rather tenuous intergalactic medium. At the same 
time, other scenarios not involving matter have also been con- 
sidered. In these cases the SMBBH is considered to be inspi- 
ralling in vacuum but in the presence of an external magnetic 



field which is anchored to the circumbinary disc (Palenzuela 
et al. 2009; Mosta et al. 2009). The extensive analysis of Mosta 
et al.[( 2009)1, in particular, has shown that, even though the elec 



tromagnetic radiation in the lowest 6 = 2 and m = 2 multipole 
accurately reflects the gravitational one, the energy emitted in 
EM waves is 13 orders of magnitude smaller than the one emit- 
ted in GW, thus making the direct detection of the two different 
signals very unlikely. The situation changes in the post-merger 
phase. In this case, in fact, the EM counterpart is supposed to be 
mainly due to the radiation from the circumbinary accretion disc, 
and, because of that, it will contain an imprint of any strong dy- 
namical change produced on the disc by the merger event. There 
are indeed two major such dynamical effects. The first one is the 
abrupt reduction of the rest-mass of the binary, emitted away in 
GWs, which is a function of the binary mass ratio, amounting 



up to - 10% for equal-mass spinning systems (Reisswig et al. 
2009) 1. The second one is the recoil velocity of the merged sys- 



tem, resulting in a kick velocity of the resulting black hole with 
respect to the hosting galaxy (see Bekenstein](| 1 973}, |Redmount| 



& Rees| ( fl989[ >) for a first discussion of the process and |Rezzolla| 
(2009|for a recent review). Leaving aside possible problems due 
to the actual value of the kicked velocity, which in some cases 
could be even larger than the escape velocity, it is clear that both 
events mentioned above can significantly affect the dynamics of 
the circumbinary disc, mainly as they contribute to the forma- 
tion and propagation of shocks, thus enhancing the possibility of 
a strong EM counterpart. 

After the first smooth-particle-hydrodynamics approach to 
the dynamical evolution of circumbinary discs performed 
by Artymowicz & Lubow (1994i, several additional numeri- 



cal investigations have been proposed in the very recent past. 
MacFadyen & Milosavljevic (2008), for example, performed 
two dimensional hydrodynamical simulations and studied in de- 
tail the evolution of the binary separation and of the disc eccen- 
tricity. By perturbing Keplerian orbits of collisionless test par- 
ticles, on the other hand, |Lippai et al.] ( |2008) > found a clear spi- 
ral shock pattern in the plane of the disc as a response to the 
kick. By performing pseudo-Newtonian numerical simulations 
of Keplerian discs O'Neill et al. ( 2009 ) have recently questioned 
the contribution of the shocks to the expected bremsstrahlung 
emissivity, while Megevand et al.| ( |2009[ ) showed that the in- 
tensity of bremsstrahlung luminosity is not much affected by 
the magnitude of the kick velocity, provided this is less than 
the smallest orbital speed of the fluid. Although they represent 
the first fully general relativistic calculations of this process, the 
simulations of Megevand et al. (2009 ) used unrealistically small 
discs which were also placed extremely close to the recoiling 
black hole. As a considerable improvement over all the previ- 
ous investigations, |Corrales et al. ( 2009 1 have carried out a sys- 
tematic study of the effects of the mass-loss and recoil over a 
number of a-discs in Newtonian gravity and two-dimensions. 
While confirming the existence of spiral shocks, they also pro- 
vided a first realistic estimate of the resulting enhanced lumi- 



nosity, which can be as large as few xl0 43 erg/s when the disc is 
assumed to be extremely efficient in radiating any local increase 
of the temperature. Very interesting results have also been ob- 
tained by Rossi et al. (2010), who estimated the maximum disc- 
to-hole mass ratio that would be stable against fragmentation due 
to self-gravity to be M c //M ~ 6 x 10~ 4 for a supermassive black 
hole with mass M = 10 6 M Q . In addition, by performing three- 
dimensional but Newtonian SPH simulations of geometrically 
thin discs, they found that the emitted luminosity corresponding 
to such small disc-to-hole mass ratios is unlikely to make the 
EM counterpart visible via wide-area sky surveys. 

In this paper we present the results of two-dimensional rel- 
ativistic numerical simulations of extended circumbinary discs 
in the post-merger phase of the merger, when the disc reacts to 
the mass loss of the central black hole and to the received kick 
velocity. By accurately capturing the dynamics of the perturbed 
disc in the relativistic regime, we investigate the dependence of 
the accretion rate on the black-hole spin and on the kick veloc- 
ity. At the same time, we introduce a new technique to locate the 
shocks that are potentially produced by the recoil and can there- 
fore assess under what conditions a spiral pattern can develop, 
producing a variability in the accretion rate and, hence, in the 
luminosity. Our "shock detector" is based on the analysis of the 
initial states of the Riemann problem solved at each cell inter- 
face and can therefore determine the location of the shock with 
extreme precision, thus revealing that the previously proposed 
criteria for the occurrence of the shock are often inaccurate. 

To compare with the general-relativistic calculations per- 
formed by Megevand et al. ( |2009 1, our initial models consider 
small-size discs with an inner radius at r ~ 40M and an outer 
one at r ~ 120M. In addition, however, we also study the dy- 
namics of large-size discs with an inner radius at r ~ 400M and 
an outer one at r ~ 4700M. These configurations have almost 
Keplerian distributions of angular momentum and are therefore 
closer to what is believed to be a realistic configuration for a cir- 
cumbinary disc. Furthermore, because the mass in the discs is 
always much smaller than the mass of the black hole (i.e. with a 
mass ratio ~ 10~ 3 ), we solve the equations of relativistic hydro- 
dynamics in the fixed spacetime of the final black hole. 

At first sight it may appear that the use of general-relativistic 
hydrodynamics is unnecessary when simulating astrophysical 
systems such as the ones considered here and especially for the 
case of large-size discs. Such a view, however, does not take into 
account that much of the dynamics in these EM counterparts 
takes place near the black-hole horizon, where general relativis- 
tic effects are not only large but essential for a correct physical 
description. Moreover, we do not have any firm theoretical basis 
to exclude small-size discs and for which the relativistic cor- 
rections are non-negligible. Finally, even in a scenario in which 
gravity could be approximated by the Newton law, we cannot 
exclude the importance of special relativistic effects. 

As all of the above mentioned investigations, also our ap- 
proach suffers from the absence of a a fully consistent treatment 
of the radiation transfer, thus allowing only for tentative con- 
clusions about the energetics involved in circumbinary accretion 
discs. However, we extend to the relativistic framework the strat- 
egy reported in Corrales et al.| ( )2009) > of performing an isother- 
mal evolution as a tool to extract luminosity curves more realistic 
than those obtained from thermal bremsstrahlung, although it ex- 
aggerates some features of the dynamics, such as the formation 
of shocks. 

The paper is organized as follows. In Section [2] we pro- 
vide the essential information about the numerical code adopted 
in the simulations. Section [3] describes the physical proper- 
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ties of the initial models, while Section |4] highlights the most 
relevant diagnostic quantities used in the rest of the paper. 
Section [5] is devoted to the presentations of the results, and, fi- 
nally, Section [6] contains a summary of our work. We assume a 
signature {-, +, +, +} for the space-time metric and we will use 
Greek letters (running from to 3) for four-dimensional space- 
time tensor components, while Latin letters (running from 1 to 
3) will be employed for three-dimensional spatial tensor compo- 
nents. Moreover, we set c — G — 1 and we extend the geometric 
units by setting m p /kB = 1, where m p is the mass of the proton, 
while kg is the Boltzmann constant. In this way the temperature 
is a dimensionless quantity. 

2. Numerical methods 

In the stationary spacetime of a Schwarzschild or Kerr black hole 
we consider the time evolution of a perfect fluid described by the 
usual energy momentum tensor 



: ph. U^U V + pg 



(1) 



where is the four velocity of the fluid, g^ Y is the space-time 
metric tensor, p is the rest-mass density, h - 1 + e + p/p the 
specific enthalpy (including rest-mass energy contribution), e the 
specific internal energy, p the thermal pressure, related to p and 
e through the usual ideal-gas equation of state (EOS) 



p=pe(y- 1) , 



(2) 



where y is the (constant) adiabatic ratio of the gas. We solve the 
corresponding equations of general relativistic non-dissipative 
hydrodynamics through the ECHO code ( Del Zanna et al.|2007) . 
Because the dynamics of the EM emission takes place on a 
timescale which is of the order of the orbital one and because 
the latter is much shorter than the viscous timescak^j] the use of 
inviscid hydrodynamics is indeed a very good approximation. 

ECHO adopts a 3 + 1 split of spacetime in which the space- 
time metric is decomposed according to 

ds 2 = -a 1 At 2 + y u (dx ; + 0dt)(dx j + f3 j dt), (3) 

where a is the lapse function, ft' is the shift vector, and y,j is 
the spatial metric tensor. The general-relativistic hydrodynami- 
cal equations are written in the following conservative form 



dfU + diT' = S, 



(4) 



which is appropriate for numerical integration via standard high- 
resolution shock-capturing (HRSC) methods developed for the 
Euler equations. The conservative variables and the correspond- 
ing fluxes in the i direction are respectively given by 





D 




aVD-ftD 












(5) 




u _ 




. aS 1 -p'U . 





whereas the sources, in any stationary background metric, can 
be written as 



5= Vr 



o 



\aW ik djy ik +Sidj0 - Udja 
iW^fidfYa + Wt J dj0 -S j dja 



(6) 



1 We recall that in geometrically thin accretion discs the local viscous 
timescale is given by t vis = r 2 /(&H 2 £l), where a is the standard alpha 
parameter ad H is the half-thickness of the disc. 



where only purely spatial quantities are present. We note that 
Vy = -^—g/a is the determinant of the spatial metric. The rela- 
tion between the evolved conservative variables (D,S j, U) and 
the primitive variables is given by 

D=pT, (7) 
Si=phT\, (8) 
U = phF 2 - p, (9) 

where T = (1 - v 2 ) -1 ^ 2 is the Lorentz factor of the bulk flow with 
respect to the Eulerian observer associated to the 3 + 1 splitting 
of the spacetime, and 



Wij=phT VjVj+pyjj, 



(10) 



is the fully spatial component of the energy-momentum tensor. 
In our setup for two dimensional disc simulations we assume the 
Kerr spacetime metric in Boyer-Lindquist coordinates (i.e. only 
ffl + 0), with the limiting case of Schwarzschild metric for van- 
ishing black-hole spins, and lay our coordinates in the equatorial 
plane of the disc (i.e. 9 = n/2). The radial numerical grid is dis- 
cretised by choosing N r points from to r max , non-uniformly 
distributed according to the following scheme 



n = r min +ai tan (a 2 x,-) 

%i ~ \Fi ~ ? min)/ (T max ^ min) 



(ii) 

(12) 



where a\ = (r max - r m i n )/ao, ai - arctanao, while r,- are the co- 
ordinate points of the uniform grid from rmin to r max . In practice, 
the free parameter ao controls the extent to which the gridpoints 
of the original uniform grid are concentrated towards r m ; n , and 
we have chosen «o = 5 in most of our simulations. The actual 
value of N r depends on the size of the disc, and it varies be- 
tween N r = 600 and = 1200. Outflow boundary conditions 
are adopted both at r m j n and r max . The azimuthal grid extends 
from to 27r, with periodic boundary conditions, and N$ = 200. 
All runs are performed with a Courant-Friedrichs-Lewy coeffi- 
cient CFL =1/2. 

The set of hydrodynamics equations is discretised in time 
with the method of lines and the evolution is performed with 
a second-order modified Euler scheme. A fifth-order finite- 
difference algorithm based on an upwind monotonicity preserv- 
ing filter is employed for spatial reconstruction of primitive vari- 
ables, whereas a two-wave HLL Riemann solver is used to en- 
sure the shock-capturing properties (see Del Zanna et al. 2007 



for further details). The timestep is generically chosen to be suf- 
ficiently small so that the second-order truncation error in time 
is comparable with the fifth-order one in space. 

As a final remark we note that as customary in HRSC meth- 
ods, we introduce a tenuous and static "atmosphere" in the re- 
gions of the fluid outside the initial model for the disc and fol- 
low the prescription detailed in Baiotti et al. (p005) for its evo- 
lution. In practice we set to zero the velocity field and reset to a 
pre-defined floor value the rest-mass density of any cell whose 
density falls below the chosen threshold value. Such a threshold 
is set to be 8 orders of magnitude below the maximum rest-mass 
density and we have checked that essentially identical results 
are obtained when changing this value of one or more orders of 
magnitude. 

3. Initial models 

As initial models we adopt stationary and axisymmetric con- 
figurations that are consistent solutions of the relativistic Euler 
equations and describe a fluid in sub-Keplerian rotation around a 
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Kerr black hole of prescribed mass and spin ( Abramowicz et al. 
1978] >. The resulting discs are geometrically thick and they can 



either have a constant or, more generally, a non-constant radial 
distribution of the specific angular momentum i. In our simula- 
tions we have considered both "small-size" models, for which 
we adopt a constant distribution of t, and "large-size''' models, 
with a distribution of I that, on the equatorial plane, obeys a 
power law 



l(r, = n/2) = Sr q , 



(13) 



where S is chosen to be positive, thus providing a disc rotation 
that is prograde with respect to the black hole rotation. A detailed 
description of the equilibrium models for non-constant specific 
angular momentum discs can be found Daigne & Font (2004). 
In particular, we have chosen a value of S such that the resulting 
thick discs possess two well defined Keplerian points, namely 
the "cusp" (which is where matter can accrete onto the black 
hole) and the "centre" (which is where the pressure has zero gra- 
dient); cf. Table[TJ 

When the exponent q in ( pj) is chosen close to 1/2, the 
rotation law tends to the Keplerian one, and the disc flattens 
towards the equatorial plane. In these circumstances the verti- 
cal structure of the disc can be essentially neglected and two- 
dimensional simulations are therefore indicative of the full three- 
dimensional dynamics. It is worth mentioning that discs with 
a rotation law given by ( 13 i have been the subject of a long- 
standing debate about whether they are subject to the so called 
"runaway instability" (|Abramowicz et al.||1983)l, which would 



lead to an exponentially rapid accretion onto the black hole ([Font 



& Daigne 2002b a, |Zanotti et al.||2003l |Daigne & Font| [2004; 
Zanotti et al.|20 05 1 Mon tero et al.|2007p . Because the onset and 



development of this instability depends on the response of the 
torus to the increased mass of the black hole, simulating this 
instability accurately requires also the evolution of the Einstein 
equations. Recent calculations of this type have been performed 
by |Montero et al.| ([2009 ) and reveal that the tori are indeed sta- 
ble irrespective of the angular momentum distribution, thus ex- 
cluding any role of the runaway instability in the dynamics of 
the discs simulated here. However, as we will further comment 
in Sect. |5.1.1| other non-axisymmetric instabilities are possible 
and have been indeed found. 

Table [T] reports the main properties of the models chosen, 
where the naming convention used allows to easily distinguish 
the small-size models (S*) from the large-size ones (L*), and 
where the number * in the name refers to the spin of the black 
hole, thus ranging between 0.00 and 0.99. As already com- 
mented in the Introduction, while the small-size models are par- 
ticularly suitable for investigating any effect of the black hole 
spin, the large-size models are those that are (astro)physically 
more realistic. The inner radius of these large-size models is typ- 
ically of a few hundreds of gravitational radii and represents the 
size of the cavity produced by the torque of the SMBBH as esti- 



mated from the expression deduced from Table 1 of Milosavljec 
|& Phinney| ( |2005] l 



' cavity 



117 

ryO.34 
a -\ > 



I ^ x 0.24 , 



M 



0.08 



[4q/(l+qY] 



2 n 0.42 



(14) 



where = a/OA, r]-\ = rj/OA (a and fi being the effective 
a-parameter of thin accretion discs and the radiative efficiency, 
respectively), Mem is the mass accretion rate in Eddington units 
and q is the mass ratio between the two coalescing black holes. 

By construction, the recoil velocities that can be studied in 
our setup are those contained in the equatorial {i.e. (r, </>)) plane 



and because it is much more advantageous to study the dynamics 
of the disc in a reference frame comoving with the black hole, we 
impose a net velocity field in addition to the equilibrium orbital 
one. In practice, at time t - Owe perform a Lorentz boost of the 
fluid velocity along the radial direction with = 0, thus mimick- 
ing a recoil velocity of the black hole in radial direction but with 
<p = 7i. We treat the imparted recoil velocity V% essentially as a 
free parameter ranging from = lOOkm/s to Vk = 3000 km/s, 
where the latter values are not realistic and serve here only to ap- 
preciate the disc dynamics under extreme conditions. We recall, 
in fact, that the recoil velocities in the orbital plane are expected 
to be < 450 km/s jKoppitz et al.||20Q7] |Herrmann e"raLl[20fJ7 
|Pollney et al.|2007) . 



In addition to the recoil, in some initial models we also con- 
sider the effects of the mass lost to gravitational waves and which 
we account by first computing the initial model in the gravita- 
tional potential of the full black hole mass, and then evolving it 
in the gravitational potential of the reduced mass. As a reference 
value we consider a decrease in the mass of ~ 3%, thus corre- 
sponding to that obtained from the typical merger of equal-mass 
spinning black holes with spins anti-aligned with the orbital an- 

Fig. 11). Higher 
itative changes in 



gular momentum (see |Reisswig et al. 2009 



values of the mass loss do not introduce qua 
the overall dynamics. 

A final comment is devoted to the EOS of the initial model, 
that we chose to be that of a polytrope p - Kp y , with y = 4/3 or 
y = 5/3. We recall that a peculiar property of these equilibrium 
models is that the ratio p/p, and therefore the temperature T and 
the sound speed c s , do not depend on the polytropic constant ^ 
which, on the other hand, determines the mass of the disc (see 
Rezzolla et al.|2003b Appendix B). As a result, the size of the 



torus is fixed, also the temperature is uniquely determined and 
cannot be rescaled further. The last column in Table [T] reports 
such a temperature at the centre of the torus, r c , as estimated 
from the ideal-gas EOS via the expression 



m p p 



(15) 



where, we recall, k B is the Boltzmann constant and m p the rest- 
mass of the proton. In geometric units and with m p /kB = 1 the 
transformation of the temperature from the dimensionless values 
to Kelvin degrees is given by 



1.088 x 10 



13 P 



K . 



(16) 



4. Methodology of the analysis 

In what follows we discuss in detail the physical quantities com- 
puted during the evolution either as representatives of the global 
evolution or of the local one. 



4.1. Global quantities 

In addition to the local Eulerian fluid variables, during the evolu- 
tion we also monitor a few global quantities that are very helpful 
for interpreting the main properties of the dynamics. These are: 
the rest-mass, the internal energy and the accretion rate at the 



2 The argument consists in proving that the function *F = k{ii + 1)0, 
where y = 1 + 1/n and p = &", does not depend on k. From this it 
follows that p/p = k® does not depend on k either. 



O. Zanotti et al.: EM counterparts of recoiling black holes: general relativistic simulations of non-Keplerian discs 



5 



Table 1. Main properties of the initial models. From left to right the columns report the name of the model, the black hole spin 
parameter a, the mass of the black hole, the disc-to-hole mass ratio, the power-law index q of the angular momentum distribution 
and the parameter S (only for models with non constant distribution of the specific angular momentum), the constant value of the 
specific angular momentum 6 (only for models with constant distribution of the specific angular momentum) the adiabatic index 
y, the inner and the outer radius of the tours, r; n and r out , the radius of the maximum rest-mass density r c , the orbital period at the 
radius of maximum rest-mass density r c , the maximum temperature T c , the orbital velocity |v^| = (v^) 1 / 2 at r; n . 



Model 


J/M 1 


M 


M d /M 




1 


S 


e 


y 


r m 


' out 


r c 






Iv1 to 






(AQ 
















(M) 


(M) 


(M) 


(K) 


(km/s) 


S.00 


0.00 


1.0 x 10 6 


2.0 x 10" 


-3 






8.0 


4/3 


40.0 


118.2 


59.8 


3.97 (h) 


5.4 x 10 lJ 


57000 


S.25 


0.25 


1.0 x 10 6 


2.1 x 10" 


-3 






8.0 


4/3 


40.0 


120.0 


60.0 


4.00 (h) 


5.6 x 10 9 


57000 


S.50 


0.50 


1.0 x 10 6 


2.2 x 10" 


-3 






8.0 


4/3 


40.0 


121.7 


60.2 


4.02 (h) 


5.7 x 10 9 


57000 


S.75 


0.75 


1.0 x 10 6 


2.1 x 10" 


-3 






8.0 


4/3 


40.0 


123.5 


60.4 


4.04 (h) 


5.8 x 10 9 


57000 


S.99 


0.99 


1.0 x 10 6 


2.1 x 10" 


-3 






8.0 


4/3 


40.0 


125.2 


60.6 


4.06 (h) 


5.9 x 10 9 


57000 


L.00 


0.00 


1.0 x 10" 


1.0 x 10- 


-4 


0.4950 


1.037 




4/3 


984.6 


403.6 


4713.7 


11.06 (d) 


8.7 x 10 b 


14970 


L.80.MA 


0.00 


1.0 x 10 6 


1.0 x 10- 


-4 


0.4950 


1.037 




5/3 


984.6 


403.6 


4713.7 


11.06(d) 


1.4 x 10 7 


14970 


L.90 


0.90 


1.0 x 10 6 


1.0 x 10- 


-4 


0.4955 


1.033 




4/3 


988.4 


400.9 


4760.0 


11.13 (d) 


7.8 x 10 6 


15030 



innermost radial point of the grid, each of which is computed as 



M disc = 2H 



J" yJyDdrdcp, 
J" e -\fyDdrd(f>, 
M(r m in) = -2H j a^Dv'c, 



Ejnt = 2H 



(17) 
(18) 
(19) 



Note that when computing the volume integral we consider the 
discs to have half thickness H which is assumed to be constant 
in radius, i.e. with H ~ c,/Q as in the standard thin disc approx- 
imation, with c s the sound speed and Q. the orbital velocity. 
In addition to ( 17 i-( 19 1 we also compute a few more diag- 



nostic quantities that, on the contrary, rely on simplified assump- 
tions reflecting the fact that the implementation does not account 
for processes such as radiation transfer and viscous dissipation. 
In particular, we compute the bremsstrahlung emissivity of the 
electron-proton collision as (Rybicki & Ligh tman|1986| l 



e BR - 2.0 x 10- 21 T l/2 Zfn e ni 
7.14 x 10 20 r 1/2 pc gs erg 



erg cm 3 s 1 



-3 -1 

cm s , 



(20) 
(21) 



where n e and «, n e are the number densities of electrons and 
ions (protons), respectively, while T is the equilibrium tempera- 
ture of both electrons and protons. The bremsstrahlung luminos- 
ity is then obtained after performing the volume integral 



L B r^ 3 x 10 78 I (T l >Yr^d*x 



Me 
M 



) erg/s 



(22) 



where the large multiplicative factor comes from the fact that 



both T and p in (22 1 are expressed in geometrized units. 



4.2. A relativistic "shock detector" 

An obvious expectation, which has been confirmed by all of the 
numerical simulations to date, is that the as the recoiling black 
hole will move in the plane of the accretion disc it will introduce 
spiral shocks which will move outwards on a timescale which 
is comparable with the orbital one. Because determining the ac- 
curate position of the shocks is important to correlate the latter 
to the EM emission, a number of suggestions have been made 
in the literature, which have a varying degree of precision. In 



particular, |Lippai et al. (2008); O'Neil l et al.| (|2009); Megevand 
|et al.| ( |20~09| ), all just looked at density and/or pressure gradients 
to infer the propagation of a spiral caustic and, therefore, of a 
possible shock (we note that in the collisionless particles treat- 
ment of Lippai et al. ( 2008 1, the existence of a shock is purely 



indicative as no shocks can be produced in this approximation). 
On the other hand, Ros si et aTT] ((2010 > used the introduction of 
an artificial viscosity, which is itself related to local density in- 
creases, to identify the location of shocks. Finally, [Corrales et al.| 
(2009) used a shock detector present in the FLASH code, which 
marks a given region as a shocked one if V ■ v < and if the pres- 
sure difference between the monitored zone and at least one of 
its neighbors exceeds the difference expected from the Rankine- 
Hugoniot jump condition for a shock of a pre-specified mini- 
mum Mach number. While more robust than those considered 
by the other authors, also this prescription is a delicate one as 
we will discuss in Sec. 15.1.11 

All of the methods mentioned above contain rather empirical 
criteria and can fail to detect shocks unless they are very strong. 
To improve the determination of the location of the shock, even 
when the latter are arbitrarily weak, we have devised a relativis- 
tic "shock detector" which exploits an idea discussed in all its 
details in |Rezzolla & Zanotti| ( |2002] > and |Rezzolla et al.| ( |2003) , 
and which consists essentially in the possibility of predicting the 
outcome of the wave pattern in a Riemann problem. (We note 
that a similar detector can be prescribed also for non-relativistic 
flows; the interested reader can find a detailed discussion in § 100 
of |Landau & Lifshitz| | |2U05) .). 

To illustrate the logic of our shock detector let us suppose 
that along a given direction, say the x-direction, two adjacent 
fluid elements 1 and 2 manifest a jump in the hydrodynamical 
quantities, such as pressure, density and velocity, thus repro- 
ducing the typical conditions of a local Riemann problem. In 
the absence of magnetic fields, the time evolution of a Riemann 
problem consists in the propagation along opposite directions of 
two nonlinear waves, either rarefactions or shocks, separated by 
a third wave, the contact discontinuity. As a result, a shock front 
will be produced if the wave pattern generated by the Riemann 
problem contains at least one shock wave, while the other wave 
can be a rarefaction wave. As shown by |Rezzolla & Zanottr] 
(2002 ), there is a simple criterion for predicting the occurrence 
of a wave pattern containing a shock wave and this amounts to 
the requirement that the relative velocity between the two states 



6 



O. Zanotti et al.: EM counterparts of recoiling black holes: general relativistic simulations of non-Keplerian discs 



1 and 2 (i.e. between two adjacent fluid cells) is larger than a 
threshold value 



Vl2 > (Vl2) s 



(23) 



where (vi2) SJ! is a function of the thermodynamic states of 1 and 
2, while V12 = (vi-V2)/(l— V1V2) is the specialrelativistic expres- 
sion for the relative velocity. When there are nonzero velocities 
in the direction tangential to the discontinuity front, the analytic 
form of (vi2) sfl is given by ( Rezzolla et al.|2003 1 



(vi2) sfi = tanh 



/ 



■ P 2 A //l 2 +^(l 



C?) 



(h 2 



- SGft> C S 



-dp 



(24) 



where J{\ = hi W\v\ while c s is the sound speed. If, on the con- 
trary, the relative velocity v l2 is smaller than (vi2) JS , then no 
shock wave can be produced and the wave pattern of the cor- 
responding Riemann problem consists of two rarefaction waves 
propagating in opposite directions. 

With this idea in mind we have constructed a sensitive shock 
detector to locate the regions of the disc which are producing a 
spiral shock. In practice we first select the direction along which 
we want to monitor the propagation of shock waves. Secondly, 
since p3) and p4| have been derived in a flat spacetime, we 
project the velocity field v ; in a local tetrad in Boyer-Lindquist 
coordinates so as to obtain the new components v J 

(25) 



r V 



V r = 



(26) 



Thirdly, we calculate v and v y from v and \r through a simple 



rotation. Finally, we compute the integral ( 24 1 in terms of the 
hatted Cartesian components and compare the result with V12. 

Note that the integral (24 1 effectively provides the minimum 
value for the occurrence of a wave pattern containing a single 
shock wave. In the limit of (vi2) s/i — > V12, in fact, the pressure 
jump across the shock wave becomes vanishingly small and a 
single rarefaction wave joining p\ and p2 propagates in the di- 
rection opposite to that of the vanishing shock wave. Therefore, 
when computing ( 24 1 we are actually integrating inside the rar- 
efaction wave, that is notoriously a self-similar solution and 



hence isentropic. This means that in evaluating (24 1 we can use 
the isentropic expression for the sound speed 



7(7 - 1)P 



\ (7 - !)P + 7P 



(27) 



where the density p is given in terms of p from p = p\(plpi) y - 

The procedure described above is completely general and 
can be proposed as an efficient shock detector for numerical rel- 
ativistic hydrodynamics. However, two subtleties should also to 
be taken into account. The first subtlety is that, for more compli- 
cated spacetimes or coordinates systems, the flat-spacetime pro- 



jection (25 1-(26 1 should be replaced by the more general form 



M'jV-i, 



(28) 



with M' j given by (Pons et al. 



1998) 



M' 








-Y [2 Yl2+y' 3 Y23 -y° ^Y22Yv-(Yn) 2 ' 



s/712 



Y22 

Vy22y33-(y23) 2 



723 
\/Y22 



\jj22 



(29) 



The second subtlety concerns the fact that, because the shock 
detector validates the inequality ( 23 1, it can be arbitrarily sen- 
sitive. Although this certainly represents an advantage, one of- 
ten wishes to disregard the whole class of weak shocks, for 
which the cont ribution to th e entropy jump is of higher order and 
As oc (A/?) 3 ( Thorne|[l973 1. In these cases the weakest shocks 
can be filtered out by making the condition < [23] > somewhat more 
restrictive and require therefore that a shock is detected if 



V12 > V12 = (vi 2 ) ss +x[(vu) ls - (V\2) SR \ 



where 



(Px-P2)(1-v x 2 V s ) 



(v s -vi){h 2P2 (r 2 m 



+ P\ -P2i 



(30) 



(31) 



with V s being the velocity of the shock wave propagating to- 
wards state 2 (see |Rezzolla et al.|2003] for the explicit expres- 
sion). Because (vi2) 2S > (vi2) ss , any value of x between and 1 
will effectively raise the threshold for the detection of the shocks, 
filtering out the weakest ones; the shocks encountered in the sim- 
ulations reported here were all rather weak and we have there- 
fore always used y — 0. The whole procedure is repeated for 
as many directions as the dimensions of the problem. Finally, 
a prescription of the relativistic shock detector as adapted for 
Newtonian fluids is presented in Appendix [A] 

5. Results 

5. 1 . Small-size models 

Although the small-size models are not astrophysically very re- 
alistic as they presume the existence of small tori in equilibrium 
near the recoiling black hole, they serve to set a comparison 
with the other general-relativistic calculations of Megevand et al. 
(2009), where similar tori were considered. In addition, by be- 



ing so close to the black hole, they are helpful in capturing those 
features of the dynamics that are most influenced by the regions 
of strong gravity. However, because of their limited extensions 
and high densities/temperatures (as an example, the model S . 08 
has p c = 3.38 x lO^g/cm 3 and T c = 7.9 x 10 & K) they will 
not be used to draw any conclusion on the emitted luminosity, 
which will be instead discussed in more detail when analyzing 
the large-size models in Sec. 5.2 



5.1.1. Shock Dynamics 

The different panels in the left column of Fig. [T] show the rest- 
mass density at three different times (i.e. t = 6.07,47.90 and 
99. 46 h) for model S. (58 and a recoiling velocity Vk = 300km/s. 
Although the imparted velocity is rather small (but close to the 
maximum possible in the orbital plane), the disc undergoes large 
variations in size and density, with a shock front that expands 
from the inner parts of the disc in an initially axisymmetric man- 
ner. This is essentially due to the reduction in the black-hole 
mass and which moves all of the equilibrium orbits to larger 
radii. As the influence area of the black hole becomes larger and 
the orbital velocities become comparable with that of the recoil, 
the disc develops shocks with the characteristic spiral structure 
discussed also in previous works (Lippai et al. 2008; Corrales 
etal.[2 009 ; Rossi et al. 2010 1 and that transports angular momen- 
tum outwards. This is shown by the panels in the right column 
of Fig. [T] which report the location of the shocks as obtained 
with the procedure illustrated in Sec. 4.2 More specifically, they 



show the quantity S d = max{0, v\ 2 - v 



}, whereby any 
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Fig. 1. Rest-mass density distributions (left columns) and shock structure (right columns) and at three different times (i.e. t = 
6.07, 47.90 and 99.46 h) for model S . 80 and a recoiling velocity V\ = 300 km/s. Note that the last panel refers to almost 25 orbital 
revolutions. The rest-mass density is plotted on logarithmic scale and in cgs units, while the shock structure is obtained by plotting 
the quantity S d (see beginning of Sec. 5.1.1 for a definition); shock waves can form in regions where S j > 0. Note that is very hard 
to locate a shock by simply looking at the density distribution. 



O. Zanotti et al.: EM counterparts of recoiling black holes: general relativistic simulations of non-Keplerian discs 



Time = 99.46 hours V_k = 300 km/s 



; ( (*f \ ; 



Time = 99,46 hours V_k = 300 km/s 





t (days) 



Fig. 2. Sound speed normalized to the kick velocity c s /Vk (left panel) and relativistic Mach number M (right panel) at t = 99.46 h 
for model S . 00 when subject to a recoil of V\ — 300 km/s. Note that Vk > c s is not a good criterion for the localization of the shock 
(cf. left panel) and that no obvious correlation is present between the supersonicity of the flow and the appearance of the shock 
(cf. right panel). 



positive value of S d marks a shocked region. Note that the re- 
gion very close to the black-hole horizon, namely at radii smaller 
than ~ IQM, is always a highly shocked one. Furthermore, as the 
evolution proceeds and the disc expands first in response to the 
mass loss and subsequently to the shocks, very little (if any) of 
the computational domain is filled by atmosphere, thus removing 
de-facto any role it can play in the dynamics of the disc. 

Figure|2]provides additional information about the dynamics 
of the disc by showing the local sound speed normalized to the 
kick velocity, and the relativistic Mach number, the latter com- 
puted as M = Fv/r s c s , where F s = 1/^1 — c\ is the Lorentz 
factor of the sound speed (Konigl 1980). Our results shows that 
the criterion suggested by Corr ales et al.| (2009) for the occur- 
rence of shocks, namely Vk > c s , can be rather misleading in the 
relativistic context. Indeed, as shown by the snapshots at time 
t - 99.46 h in Fig. [T] and which refers to almost 25 orbital rev- 
olutions, a clear spiral shock forms even if the sound speed is 
more than 30 times larger than the kick velocity. Indeed, the left 
panel of Fig.[2]seems to suggest that if any, c s /Vk » 1 is possibly 
a reasonable necessary condition for the approximate location of 
the shock. In addition, the correlation between the occurrence of 
a shock and the local sound speed is very weak and this is appar- 
ent in the right panel of Fig. [2] where it is clear that the flow is 
highly supersonic in the inner regions of the disc and mildly sub- 
sonic in the outer regions. Yet, the spiral-shock structure extends 
continuously across the whole disc. We also note that although 
the precise morphology of the spiral shocks will depend on the 
spin of the black hole, this dependence is only very weak and all 
the considerations made above for model S . 00 hold true quali- 
tatively also for spinning black holes. 

It is finally worth remarking that the shocks formed here 
are very mild and not relativistic. Even for Vk = 3000 km/s, 
the shock velocity maintains a typical value V s ~ 0.15 and the 
velocity jump at the shocks is also rather limited, producing a 
v/Av ~ 30. This means, for instance, that such shocks are un- 
able of accelerating electrons through the classical mechanism 
of |BelH < |l978) . 




0.04 0.06 0.08 

frequency (mHz) 



0. 1 



Fig. 3. Time evolution of the internal energy when normalized 
to the its initial value (top panel) and the corresponding power 
spectrum (bottom panel) in a model with V\ = and with a mass 
loss of 1% the initial mass of the black hole. 



5.1.2. Mass loss and Quasi-Periodic Dynamics 

The dynamics of the disc can change considerably if the black 
hole is assumed to be recoiling with negligible velocity in the 
orbital plane and only mass loss is taken into account. By con- 
sidering mass losses in the range 1%-10%, |Q'Neill et al.| ([2009) 
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Fig. 4. Mass accretion rate measured at r — r m ; n for different Fig. 5. Mass accretion rate measured at r — r m j n for the small- 
values of the black hole spin-parameter in the small-size models, size model S . 0® and different values of the recoil velocity. 
i.e. S . 0Q-S . 99. Shown in the inset as a function of the black- Shown in the inset as a function of the recoil velocity is the rel- 



hole spin is the stationary accretion rate reached after 5 d. 



ative baryon mass accreted after 3 d. 



showed that shocks can form even in the absence of a recoil ve- 
locity, provided that the mass loss is larger than the half thick- 
ness of the disc. The perturbation induced by the mass loss is 
spherically symmetric and it causes the disc to expand as each 
fluid element will want to move to the larger radii corresponding 
to the equilibrium orbit for the given initial angular momentum. 
Together with this expansion, however, restoring forces will also 
induce the disc to contract in the effective-potential well of the 
black hole with the characteristic frequency of the lowest or- 
der p-mode, and which is not too different from the epicyclic 
frequency at the disc centre (Rezzolla et al. 2003b). We recall 
that the restoring force responsible for the appearance of such 
p modes is a combination of pressure gradients, centrifugal and 
gravitational forces, with the last two playing the dominant roles 
for the discs considered here (Kato 2001 ; Rezzol faet al.|20 03b). 

The oscillating behavior induced by the sudden change of 
the potential well and the subsequent development of the insta- 
bility is shown in Fig. [3] where the top panel reports the time 
evolution of a typical global quantity, i.e. the internal energy, 
when normalized to its initial value. Interestingly, the remark- 
able periodicity that characterizes the dynamics is the same as 



found by Zanotti et al. ( 2003| > when studying global modes of 
oscillation of thick discs around black holes. The corresponding 
power spectrum is shown in the bottom panel of Fig. [3] obtained 
through a FFT of the time series for t < 3 d reveals the presence 
of a fundamental mode of oscillations at / ~ 4.17 x 10~ 5 Hz and 
of two overtones. The first overtone is at o\ ~ 6.28 x 10~ 5 Hz, 
while the second one, very close to twice the fundamental fre- 
quency, 02 ~ 8.37 x 10~ 5 Hz ~ 2f, is produced by nonlinear 
coupling of the fundamental mode with itself (see Zanotti et al. 



observed in the QPOs of low-mass X-ray binaries containing a 
black hole (see Remill ard & McClintock|2006| for a recent re- 
view) and for which a simple model based on basic p-mode os- 
cillations of a small accretion torus orbiting close to the black 
hole was recently proposed (Rezzolla et al. 2003a), and which 
remains one of the most convincing explanations of the observed 
phenomenology (Schnittman & Rezzolla 2006}. 

It is difficult not to note that the harmonic behaviour shown 
in the top panel of Fig. [3] is lost at t 3d and that the inter- 
nal energy increases monotonically after that. This is due to the 
onset of a non-axisymmetric instability which produces spiral 
arms that rapidly spread to cover the whole disc. An instabil- 
ity of this type was not pointed out by Megeva nd et al.| ([2009 ) 
although they have used similar models and we believe that 
this is probably because their simulations were interrupted af- 
ter ~ 1 1 h (or ~ 6 orbital periods as measured at the point of 
maximum rest mass density), which is too early for the devel- 
opment of the instability. On the other hand, such type of insta- 
bilities in non-Keplerian discs have been discussed by a number 
of authors, starting from the pioneering work by Papaloizou & 
Pringle ( 1984 ). A detailed comparison between the linear pertur- 
bative analysis of these instabilities and two-dimensional numer- 
ical simulations in a Schwarzschild spacetime was already pro- 
posed more than twenty years ago by Blaes & Hawley| ([T988 ), 
who found the development of the same spiral structures, which 
transport both mass and angular momentum outwards even in 
the absence of mass losCD While we cannot concentrate here 
on a detailed discussion of these instabilities, it is sufficient to 



2005 ). Collectively, these modes of oscillations provide a series 
of modes in the same ratio of the integer numbers 2:3:4 



3 We have verified that the spiral arms do indeed develop also in the 
absence of a mass loss or of a recoil and that also in this case the insta- 
bility takes place after ~ 4 d. 



10 



O. Zanotti et al.: EM counterparts of recoiling black holes: general relativistic simulations of non-Keplerian discs 



remark that a spiral-shock pattern and all of the associated phe- 
nomenology, can be generated even when the recoil velocity is 
zero. 

Overall, the dynamics observed for model S . 00 suggests that 
transient oscillating phenomena may exist in the post-merger 
phase of SMBBH. In this case, the occurrence of QPOs in the 
accretion and thus in the luminosity of potential EM counter- 
parts, followed then by the development of non-axisymmetric 
instabilities would be a unique and convincing signature that a 
SMBBH merger with small recoil velocities has taken place. 



5.1.3. Accretion rates 



As originally pointed out by Kozlowski et al. (1978 1, in non- 
Keplerian discs no viscosity is needed in order to support accre- 
tion in the vicinity of the cusp. Figure |4]reports the baryon mass 
accretion rate measured at r — r m j n for different values of the 
black hole spin-parameter in the small-size models, i.e. S.00- 
S . 99, and provides further support to the interpretation of the 
shock dynamics discussed above. In particular, it is easy to real- 
ize that because of the sudden mass loss and hence of the reduced 
gravitational attraction, the disc reacts to the excess of angular 
momentum by expanding. This effect only lasts for a couple of 
orbital periods after the merger and leads to the large decrease in 
mass accretion rate shown in Fig. |4]for t < 0.6 d. 

Subsequently, as the effect of the kick velocity extends to re- 
gions of the flow with smaller orbital velocities and becomes 
dominant, non-axisymmetric density structures form and the 
perturbed disc starts filling the low density central cavity while 
increasing the accretion rate. This is reflected by the the large 
increase in the mass accretion rate shown in Fig.|4]for t > 0.6 d, 
which is essentially independent of the black-hole spin. 

Since our treatment does not include the compensating ef- 
fect of the radiation drag exerted by the photons on the ac- 
creting matter, the accretion rate increases undisturbed reaching 
values that are ~ 6 orders of magnitude above the Eddington 
limit. After about 6 orbital periods (« Id) M saturates and after 
~ 5 d all of the models have lost ~ 32% of their mass. Fig. [4] 
also shows that the spin of the black hole has little influence 
on the dynamics of the disc, although not as little as inferred 
from Fig. 13 of Megevand et al. ( 2009| l. In particular, we have 
found that after ~ 3 d the accretion rate slightly decreases when 
increasing the spin-parameter (see the inset of Fig. [4}, so that 
M| a= o ~ 2.5 M\ a =o .99 after ~ 5 d of evolution. Interestingly, this 
dependence of the accretion rate on the spin of the black hole is 
rather generic and has been found also in other simulations (not 
reported here) where no perturbation on the flow is introduced. 

Of course, the effect of an increasingly large kick velocity 
on the disc dynamics is much more pronounced in this case, as 
it emerges from Fig. [5] which reports the accretion rate for dif- 
ferent values of Vk. As already found by Megevand et al. (2009 ), 
larger recoil velocities tend to anticipate the occurrence of the 
burst in the accretion rate and, simultaneously, produce larger 
absolute values of M at the burst. After the initial burst, however, 
and no later than ~ 2 days, the accretion rates become nearly 
independent of the recoil velocity. The differences in M are il- 
lustrated in the inset of Fig. [5] which shows the relative baryon 
mass accreted after 3 d and which shows a variation of ~ 50% at 
most. It is also interesting to note that there is no clear imprint 
of the spiral shock dynamics onto the accretion rate. This is due 
to the fact that the spiral shocks essentially redistribute angular 
momentum and thus modify the disc structure and dynamics far 
from the black hole. 



5.2. Large-size model 

We next switch to discussing the dynamics of the large-size 
model L.00 for different values of the recoil velocity. We re- 
call that there are at least two different reasons to consider this 
second class of models. The first one is that they reflect the ex- 
pectations of a circumbinary disc: they are quasi-Keplerian, ex- 
tended and at a large distance from the binary. The second one 
is that by having a lower density, i.e. p c = 1.38 x 10~ I0 g/cm 3 , 
and hence a lower temperature, i.e. T c = 8.7 x 10 6 K, they lead 
to much more reasonable values of the recoil-enhanced luminos- 
ity. In what follows we briefly review the overall dynamics and 
then concentrate on how to compute a realistic estimate for the 
luminosity. 

5.2.1. Shock Dynamics 

The overall dynamics of the large-size models is qualitatively 
similar to that of the smaller counterparts but with three impor- 
tant differences. The first one is a much weaker dependence of 
the dynamics on the recoil velocity. This is essentially due to the 
fact that the inner edge of the discs is so far from the recoiling 
black hole that only extremely large (and unrealistic) values of 
the recoil velocity induce a significant modification of the orbital 
velocity. This is shown in Fig. [6| which reports in the left col- 
umn the rest-mass density after 60 and 180 d for a recoil velocity 
Vk = 500km/s applied to model L . 00. The right column shows 
instead the corresponding shock structure with the same notation 
used in Fig. [T] and highlights that a clear spiral structure is lost 
already after 180d, which now corresponds to almost 18 orbital 
periods. Note that the temperature distribution is inversionally 
proportional to that of the density (not shown in Fig. [6} and that 
the central region is filled only with the atmosphere fluid and 
thus appears as white in the left column. However, several small 
shocks are produced in this cavity and these are clearly revealed 
by the shock detector images on the right column which reports 
the central region as dark; this is just an artifact that does not 
have a dynamical impact. The behaviour in Fig.|6]should also be 
contrasted with the spiral-shock structure of model S . 00, which 
instead persisted intact after almost 25 orbital periods (cf. bot- 
tom right panel of Fig. [T}. This behaviour and the rapid disap- 
pearance of the spiral shock structure is further pronounced as 
the recoil velocity is increased to Vk = 1000 km/s (not shown 
here). 

The second difference is that for sufficiently large recoils, 
namely for Vk > 2000 km/s, the initial cavity between the black 
hole horizon and the inner edge of the disc can be filled rapidly 
by infalling material. This is evident when looking at Fig. [7] 
which reports the rest-mass density and shock structure after 60 
and 1 80 d for a recoil velocity Vk = 3000 km/ s applied to model 
L . 00. When contrasted with Fig. [6] which has the same spatial 
extent, it is easy to notice that already after ~ 60 d, or ~ 5 orbital 
revolutions, the central cavity is filled with high-density mate- 
rial, some of which extends right onto the black hole. Similarly 
to what already seen for smaller recoils, also the right column of 
Fig. [7] shows that in this case the spiral structure is not present, 
although spatially extended shocks are formed both in the in- 
ner regions and in the outer parts of the disc. Note that the ve- 
locity jump at the shocks is again rather small, being at most 
Av ~ 2.5x 10~ 4 , hence insufficient to accelerate particles through 
the various acceleration mechanisms involving shock waves. 

Interestingly, because the black hole is moving at very large 
velocities in an ambient fluid which also has a non-negligible an- 
gular momentum, a low-density region which resembles a horn- 
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Fig. 6. Left column: Rest-mass density after 60 and 180 d for a recoil velocity V± = 500 km/ s applied to model L . 00. The scale is 
logarithmic and expressed in cgs units. Right column: Shock structure as presented in Fig.[T]for the same panels in the left column. 
Once again we remark that is very hard to locate a shock by simply looking at the density distribution, especially when they are 
very weak (here we have used^- = in Eq. |30|l. Note also that the temperature distribution is inversionally proportional to that of 
the density (not shown here). 



shaped "cavity" is produced in the downstream part of the flow 
when Vk = 3000 km/s at t — 180d after the merger. This cavity 
is shown in Fig. [8] and its orientation is directly related to the 
direction of the recoil velocity. A simple change of sign in the 
recoil velocity, in fact, would rotate the cavity of 180 degrees 
around the black hole (not shown in Fig. [8]). The formation of 



such a cavity is noticeable also in the simulations of |Rossi et al 
(2010), where however it is not discussed. The cavity leads to 



quasi-periodic variation of the accretion rate as clumps of matter 
in the downstream of the flow enter the cavity and streams onto 
the black hole (cf. the oscillations of M in Fig. [9] after t ~ 200 d 
for Vk = 3000 km/s). The generation of this flow pattern has an 
interest of its own, being a non-trivial variant of the Bondi-Hoyle 
accretion flow onto a moving black hole and will be investigated 
with greater detail in a distinct work. 

The third and final difference is really a combination of the 
phenomenology discussed above as it manifests itself in terms 
of the accretion rate. This is shown in Fig. [9] which reports the 



mass accretion rate at r — r m j n for different values of the kick 
velocity in the large-size model L . 00. As already discussed with 
Fig.[5]for the small-size models, also in these larger discs the ac- 
cretion rate shows a rapid increase when the black hole "meets" 
the disc, reaching values M 10M Q /yr, that are well above 
the Eddington one because of the absence of the radiation-drag 
contribution. Differently from the smaller discs, however, the lag 
between the merger and the increase of the accretion rate is not 
linear but rather increases nonlinearly as the recoil velocity is 
decreased. This is clearly shown in Fig. [9] where the accretion 
rate jumps to high values after ~ 35, d for Vk = 3000 km/s, 
after ~ 100, d for Vk = 2000 km/s and after almost one year, 
i.e. ~ 330, d for Vk = 1000 km/s. The reason for this nonlinear 
response of the disc has to be found in the fact that for smaller re- 
coil velocities the time required by the black hole to cross the ini- 
tial cavity T cross will be much larger than the typical orbital rev- 
olution time. As an example, for Vk = 1000 km/s the timescale 
ratio is t cioss /t c ~ 33, so that the disc will have sufficient time to 
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Fig. 7. The same as in Fig.^but for a recoil velocity Vk = 3000 km/s. Note that the spiral-shock structure is never present and that 
the inner cavity is rapidly filled by accreting gas. In addition, an oblique shock comprising a low-density region is formed in the 
inner parts of the flow. 



readjust itself to the new gravitational field and thus redistribute 
its orbital angular momentum. In practice the inner edge of the 
disc will move to larger radii as a result of the mass loss and of 
the varied potential, thus delaying nonlinearly the contact with 
the black hole and thus the steep increase in the accretion rate. 

As a final remark we note that all of the phenomenology dis- 
cussed here for the model L . 06) has been found also for a large- 
size model around a rapidly spinning black hole, namely L . 90. 
Because the overall differences are minute and of the order of 
few percent at most, they will not be discussed here. The reasons 
behind these similarities are rather obvious: the large-size discs 
have inner radii that are too far from the black hole to be sen- 
sitive to the spin-induced corrections which decay much more 
rapidly and as 1/r 3 . 



5.2.2. EM Luminosities 



thought to be an efficient emission mechanism through which 
circumbinary discs may become visible in the electromagnetic 
spectrum (|Megevand et al.|2009||Corrales et al.|2009||Bode et al. 
2009; Anderso n et al.|2009[ ). However, fhermal-bremsstrahlung 
emission from circumbinary is affected by a serious problem 
which has been so far underestimated or not sufficiently em- 
phasized. This has to do with the fact that bremsstrahlung cool- 
ing time is too short (Corrale s et al.|200 9) or, stated differently, 
that the internal energy budget of the emitting gas is not large 
enough to allow for the bremsstrahlung emission to last but for a 
few seconds. This can be easily estimated as f coo i = E^/Lbr, 



with Ei nt and Lbr obtained from (18 1 and (22 1, respectively. 
For the large model L.00 we have E ml ~ 3.4 x 10 50 erg and 
Lbr ~ 2.8 x 10 49 erg/s at time t — and this estimate re- 
mains of the same order of magnitude during the evolution. This 
yields to f coo i 12 s. The situation is even worse if we consider 
the transition to the relativistic regime. In this case, in fact, not 



Because of the relatively high temperature of the gas and of 
the generation of a shock pattern, thermal bremsstrahlung is 
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Fig. 8. Magnification of the central region of the rest-mass den- 
sity in model L.00 after 180d and for a recoil velocity = 
3000km/s. Note that the colormap is slightly different from the 
one in Fig. |7]to highlight the presence of a cavity. 
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Fig. 9. Mass accretion rate at r — r m ; n for different values of the 
kick velocity in the large-size model L . 00. 

only the bremsstra hlung emissivity is increase d by a factot^joc 
[1 +4AT/(1Q 10 K)] ( |Rybicki & Lightman|l986| , but also the col- 
lisions between particles of the same species start contributing 



4 It should be remarked, however, that when the electron become 
relativistic, i.e. for T > 5.9 x 10 9 K, other emission mechanisms, such 
as inverse Compton or synchrotron (if a magnetic field is also present), 
are generally more efficient. 



Fig. 10. Top panel: Luminosity computed in the isothermal evo- 
lution L; sot for different values of the kick velocity for model 
L.00 and for a poly tropic index y = 4/3. Note the presence 
of a peak at about ~ 20 d after the merger, of a binary with 
total mass M 10 6 M Q and of a persistent luminosity for sev- 
eral days at values which are a factor of a few smaller. Bottom 
panel: Comparison of L; sot for model L . 00 when computed with 
a polytropic index y = 4/3 (thick lines) or y = 5/3 (thin lines). 
The comparison is made for two reference recoil velocities and 
shows that the results are very similar, although a stiffer EOS 
leads to slightly larger luminosities. 



significantly to the bremsstrahlung emission (Svensson 1982) 
through radiation in moments other than the electric dipole 
(which is strictly zero for particles of the same species ( Krolik 
|T999l ». 

Of course, there are also other factors that can work in favour 
of a bremsstrahlung emission and which we have not taken 
into account. A first one is that we have neglected the thermal 
bremsstrahlung absorption, which is likely to enhance signifi- 
cantly the bremsstrahlung cooling time by acting as a source of 
additional internal energy. Moreover, it is also possible that the 
spiral shock originating from the very central region dissipates 
considerably as it propagates outwards, hence confining the bulk 
of the bremsstrahlung luminosity from within a very small por- 
tion of the disc (we recall that the bremsstrahlung luminosity 
is proportional to the volume integral over the emitting source). 
Overall, while we cannot rule out thermal bremsstrahlung as an 
emission mechanism from the circumbinary disc, it is also ev- 
ident to us that the luminosity estimates made so far without a 
proper treatment of radiation transfer are excessively optimistic. 

A second possible estimate of the luminosity is given by the 
accretion-powered luminosity, L acc = rjMc 2 , where r\ is the ra- 
diative efficiency. However, lacking any treatment of the radia- 
tion pressure effects, such an estimate would only provide mis- 
leading conclusions and therefore we will note use it hereafter. 
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A third and possibly more accurate estimate of the luminos- 
ity can be made by assuming that all the changes in the temper- 
ature that are due to a local compression will be dissipated as 
radiation. This idea, proposed in Newtonian physics by Corrales 
et al. ( 2009| l, can be summarized and implemented in a general 
relativistic context as detailed below. 

Consider the evolution of the disc with an equation of state 
p(T) = pkbT/nip and a specific internal energy given by e(T) = 
ki,T/[(y - l)m p ] = Iplp, where the last equality has been ob- 
tained for y = 5/3. In general, there is no necessity to evolve 
the energy equation in an isothermal evolution since the energy 
can be computed directly from the temperature and the latter is 
constant by construction. However, the internal energy can be 
nevertheless evolved in time with the only aim of computing the 
difference p[e - e(T)], which is then assumed to be radiated in- 
stantaneously. The relativistic equation for the evolution of the 
total internal energy density e = p( 1 + e) is ( |Anile|1990] > 



uTV^e + (e + p)& = 0, 



(32) 



where = V^u? is the expansion of the fluid. The continuity 
equation V ^(pu^) = can then be used to rewrite Eq. (32 1 as 



d,( yjyWpe) + <9 ; [ y/ypeWiav* - /?')] = 

-pd,( y/yW) - pdi(a yfyu') 



(33) 



One aspect to note is that Eq. p3f is not written in a conser- 
vative form because of the derivatives on the right hand side act- 
ing on the flow variables. While this is not ideal within our for- 
mulation of the hydrodynamics equations, the modifications are 
minimal. Indeed, since the Lorentz factor does not change signif- 
icantly during the evolution, it is reasonable to neglect the term 
oc <9f( y/yW), while the spatial derivatives of the term d,(a y/yu 1 ) 
can be treated with standard finite-difference methods without a 
significant loss of accuracy across the discontinuities. 

In practice we have assumed an initial temperature for the 
disc which is uniform in space and set it to be Tq = lT c , 
where T c is the maximum temperature at the center of the disc 
(cf. Table[T|. An estimate of the luminosity is then trivially com- 
puted by performing at each timestep a volume integration of 
the difference p[e - e(T)] and by dividing it by the simula- 
tion timestep. Finally, we reset the specific internal energy to 
e = e(To), so as to guarantee that the evolution is effectively 
isothermal. 



The top panel of Fig. 10 shows the luminosity computed in 
this way for model L . 00 with a polytropic index y = 4/3 and for 
different values of the recoil velocity. Following |Corrales et al.| 
(2009), the values reported have been computed by neglecting 
the negative contributions to the luminosity that are produced in 
regions experiencing rarefactions. Because when accounted for 
these negative contributions typically yield values that are one 
order of magnitude smaller, the values in Fig. 10 should be taken 
as upper limits to the emitted luminosity (see |Corrales et al.| 
2009| l. Clearly, the evolution of the emitted energy has a peak 
that is larger for stronger recoils and that appears at t ~ 33, 18, 
and 7 d after the merger for V k = 300, 1000 and 3000 km/s, 
respectively. While the peaks are the consequence of the strong 
shocks that are produced in the inner parts of the disc as the latter 
approaches the black hole, the asymptotic values of the isother- 
mal luminosity are instead produced by the local compressions 
in the disc. As such, the peaks in the luminosity are not related 
to the encounter of the black hole with the disc and therefore 
they are not correlated with the increase in the mass accretion 
rate, which in general takes place at later times (cf. right panel of 
Fig. [TT) . The bottom panel of Fig. [10] on the other hand, shows 



a comparison of L; sot for model L.00 when computed with a 
polytropic index y - 4/3 (thick lines) or y = 5/3 (thin lines). 
The comparison is made for two reference recoil velocities of 
Vk = 300 km/s and Vk = 1000 km/s and it shows that the results 
are very similar, although a stiffer EOS leads to slightly larger 
luminosities. 

Overall, our estimates of the luminosity computed within the 
isothermal evolution approximation confirm those by Corrales 
et al. ( 2009 1 even though the temperature of our large size model 
is one order of magnitude larger than that reported in Fig. 7 
of |Corrales et al.| ( |20 09). While we believe that these estimates 
of the luminosity are the most reasonable ones that can be ob- 
tained with a code that is intrinsically unable to account for ra- 
diation losses, we should also stress that the isothermal evolu- 
tion by itself provides a less reliable description of the overall 
dynamics. This is evident, for example, from Fig. 11 whose 
left panel offers a comparison of the rest-mass density profiles 
along <p — as obtained with the standard evolution of the en- 
ergy equation (red solid line) and with the isothermal evolution 
(blue dashed line) for the large-size model L . 00 and a recoil of 
Vk — 1000 km/s. Note that the isothermal evolution tends to in- 
crease the density gradients, especially during the initial phases, 
which is also when the peak in the luminosity appears in the 
light curves. Also shown in the right panel of Fig.[TT]is the mass 
accretion rate at r — r m ; n for the isothermal evolution and dif- 
ferent values of the kick velocity. This panel shows that since 
the discs are more sensitives to compressions, their response to 
the recoil is different and, in particular, takes place earlier than 
in the full-evolution case (cf. with Pig. [9j. Furthermore, with the 
exception of the very large kick, the mass accretion rate does not 
show a very large jump, but rather a first small jump followed 
by smooth increase to the final asymptotic behaviour which is 
reached at times comparable to those of the full evolution. This 
is clearly the case for Vk = 1000 km/s and is due to the fact that 
accretion starts earlier for this matter whose energy has been de- 
creased by the radiative losses. 

In conclusion, and in spite of the caveats made above, we 
believe that luminosities as large as few L 10 43 erg/s should 
be expected at about ~ 30 d after the merger of a binary with total 
mass M 10 6 M o and that these luminosities should persist for 
several days to values which are a factor of a few smaller. 

6. Conclusions 

We have presented the results of two-dimensional general rela- 
tivistic numerical simulations of small and extended circumbi- 
nary discs in the post-merger phase of the merger, when the disc 
reacts to the mass loss of the central black hole and to its re- 
coil velocity. Our analysis benefitted from being able to capture 
accurately the dynamics of the perturbed disc in the relativis- 
tic regime, thus allowing us to investigate the dependence of the 
accretion rate on the black-hole spin and on the kick velocity. 
Furthermore, by considering discs that are quasi-Keplerian, ex- 
tended and at a large distance from the binary, we were able to 
consider realistic scenarios even in the general relativistic frame- 
work. 

Another important aspect of our work is the use of a novel 
and accurate technique to construct a "shock detector" and hence 
to determine where, within the flow, the shocks produced by the 
recoil are located. This, in turn, has allowed us to assess un- 
der what conditions a spiral-shock pattern can develop, produce 
a variability in the accretion rate and, hence, in the luminos- 
ity. Our relativistic shock detector (for which we also present 
a Newtonian equivalent in Appendix [A]) is based on the analysis 
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Fig. 11. Left panel: Comparison of the rest-mass density along the <p — direction as obtained with the standard evolution of the 
energy equation (red solid line) and with the isothermal evolution (blue dashed line), at two different times as shown in the two 
panels. The data refers to model L.88 with a recoil of 1000 km/s. Right panel: Mass accretion rate at r = r m \ n for the isothermal 
evolution and different values of the kick velocity in the large-size model L.00. This panel should be compared with Fig. [9] and 
which refers to an evolution of the energy equation. 



of the initial states of the Riemann problem solved at each cell 
interface and can therefore determine the location of the shock 
with the same resolution as that of the spatial grid, revealing that 
the previously proposed criteria for the occurrence of the shock 
are often inaccurate. 

Overall, we can confirm within a general relativistic regime 
many of the results found previously in Newtonian or pseudo- 
Newtonian gravity. More specifically, we find that for discs 
that are sufficiently small and close to the black-hole, a regu- 
lar spiral-shock develops as a result of the recoiling black hole. 
The strength, shape and persistence of the shocks, however, de- 
pend sensitively on both the size of the tori and on the intensity 
of the recoil. As a result, while the spiral shock is stable over 
many orbital periods in the case of small discs subject to small 
recoils, it never develops or is rapidly destroyed in discs that are 
large and subject to large recoil velocities. It is worth noting that 
the typical velocity jumps at the shocks are Av < 2.5 x 10~ 4 , even 
with a kick velocity = 3000 km/s. It is therefore possible that 
such shocks may be damped by dissipative viscous processes or 
radiative losses. 

Besides the interesting shock properties described above, we 
also found that the disc dynamics is only very weakly dependent 
on the black-hole spin. The latter influences only the small-size 
tori by modifying the accretion rate and by leading to smaller 
accreted masses per unit time for more rapidly spinning black 
holes. Finally, we found that even in the limit of a vanishing 
kick velocity and as long as a mass loss is present, the disc 
goes through a phase of regular oscillations characterized by 
the excitation of the p modes of the disc; this is then followed 
by the appearance of a spiral-shock pattern generated by non- 
axisymmetric instabilities affecting the disc. This opens the door 



to the possibility that quasi-periodic oscillations are observed in 
the post-merger phase of SMBBH. 

Computing the EM counterpart to the merger event repre- 
sents of course a fundamental aspect of our investigation. In 
contrast with other works, however, we have questioned the es- 
timates of the bremsstrahlung luminosity when computed with- 
out properly taking into account the radiation transfer. The en- 
ergetic reservoir available for the EM emission is, in fact, too 
small to yield but unrealistically short cooling times. While we 
cannot rule out thermal bremsstrahlung as the main EM emission 
mechanism, we believe that the bremsstrahlung-luminosity esti- 
mates made so far without a proper treatment of radiation trans- 
fer are excessively optimistic. A somewhat realistic estimates 
of the emitted luminosity can be obtained by assuming that the 
internal-energy enhancements due to local compressions in the 
perturbed disc are immediately radiated away, i.e. by consider- 
ing an "isothermal" evolution like the one recently investigated 
in Newtonian physics by Corrales et al.| ( 2009) . In this case, and 



despite the fact that the isothermal evolution tends enhance the 
compressibility of the fluid, we find that the energy emitted can 
reach a peak value above L =s 10 43 erg/s at about ~ 30 d after 
the merger of a binary with total mass M =t 10 6 M o and persist 
for several days at values which are a factor of a few smaller. If 
confirmed by more sophisticated calculations such a signal could 
represent an interesting EM counterpart of the merger of binary 
black-hole system. 

As a final remark we note that while a rather robust picture 
is emerging from the collective work done so far on the post- 
merger dynamics of the circumbinary disc around a SMBBH, 
much remains to be done to compute realistically the resulting 
EM emission. Important improvements to the treatment consid- 
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ered here must include the presence of a magnetic field, a proper 
treatment of the radiation transfer and the extension to three- 
dimensional calculations. All of these will be the focus of our 
future work on this subject. 
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Appendix A: Newtonian version of the shock 
detector 

In this Appendix we provide the basic expressions that allow 
to build the Newtonian version of the relativistic shock detector 



described in Sec. 4.2 The logic, of course, is exactly the same 
and it consists of the following steps. 



L. 



First choose the direction, say the x-direction, along which 
to monitor the generation of shock waves. 
Given any two adjacent cells, label with 1 the cell with higher 
pressure and with 2 the other one. 

Compute the relative velocity vi 2 = V\ — V 2 along the 
^-direction and compare it with the value (see Landau & 
|Lifshitz|2004| §100) 



(Vl2)„ 



Newt 



y- 



c s (pi) 



l - 



(r-D/2y 



(A.l) 



where y is the adiabatic index of the gas and c s (p\) is the 
sound speed in the state 1. It is interesting to note that in 
the Newtonian case the tangential velocities do not affect 
the actual value of the threshold (vi 2 ) M , which, therefore, 
maintains the same expression as for one-dimensional flows. 
Indeed, the fact that in the relativistic regime the threshold 
given by ( 24 1 does depend on the tangential velocities is at 
the origin of the appearance of new relativistic effects dis- 
cussed in |Rezzolla & Zanotti| (j2002). 
A shock is therefore detected if 



Vl2 > (Vt 2 ) ss 



(A.2) 



The procedure is repeated for as many directions as the di- 
mensions of the problem. As commented in the main text, we 
note here too that it may be necessary at times to filter-out the 
smallest shocks and therefore make the condition for shock de- 
tection more stringent. In analogy with Eq. ( |30"| ), the following 
condition can then be implemented 



V 12 > V12 = (Vl 2 ) M 


+ X 




- (V12) M 




, (A.3) 




Newt 




Newt 


Newt- 





where (see |Landau & Lifshitz|2004l §100) 

(vi 2 ) 2 



Oi - pi) 



\P2 [(y+ +(y-i)/>2]' 



(A.4) 



and where the parameter x e [0, 1 ] needs to be tuned to the 
desired level of accuracy. 
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